Introduction into CAGE data

CAGE Analysis with CAGEr

Loading in data from an R-package

Here, we will do a sample workflow for two samples at different stages in zebrafish development: 512 cell stage (~2,7 hours hpf) and Prim 6 (~24 hpf). These stages have been previously published by Nepal et al (2013) and are available in an R package (by Haberle). Similar packages are available for mice and human data.
For the first part of the workflow we will use the R package called “CAGEr”. We will do some of the standard workflow and gradually manipulate the data ourselves in R.
So let’s start by loading in the package for the developmental stages as well as CAGEr:

# packages
require(ZebrafishDevelopmentalCAGE)
require(CAGEr)

Loading the samples of which we have data and display which stages we have here:

# loading data
data(ZebrafishSamples)
head(ZebrafishSamples) 

# samples
samples <- as.character(ZebrafishSamples$sample) # to see all the samples
samples

 [1] "zf_unfertilized_egg" "zf_fertilized_egg"   "zf_64cells"         
 [4] "zf_512cells"         "zf_high"             "zf_oblong"          
 [7] "zf_sphere_dome"      "zf_30perc_dome"      "zf_shield"          
[10] "zf_14somites"        "zf_prim6"            "zf_prim20"          

As said earlier, we will pick just two developmental stages: “zf_512cells” and “zf_prim6”. In the vignette (and manual) the way of loading in these data is by the function importPublicData(). We need to specify the source, dataset, group, and samples. The source is “ZebrafishDevelopment” and the rest are found in the dataframe of the ZebrafishSamples: ZebrafishCAGE, development, samplenames. The avoid spelling everything out we’ll use the vector of samples:

# creating a CAGEset object:
myCAGEset <- importPublicData(source = "ZebrafishDevelopment", dataset = "ZebrafishCAGE", 
group = "development", sample = samples[c(4,11)]  )
# what does it look like?
myCAGEset

S4 Object of class CAGEset

=======================================
Input data information
=======================================
Reference genome (organism): BSgenome.Drerio.UCSC.danRer7
Input file type: ZebrafishDevelopment
Input file names: ZebrafishDevelopment__zf_512cells, ZebrafishDevelopment__zf_prim6
Sample labels: zf_512cells, zf_prim6
=======================================
CTSS information
=======================================
CTSS chromosome: chr1, chr1, chr1, ...
CTSS position: 3783, 3788, 3798, ...
CTSS strand: +, +, +, ...
Tag count:
    -> zf_512cells: 0, 0, 0, ...
    -> zf_prim6: 1, 1, 1, ...
Normalized tpm:
=======================================
Tag cluster (TC) information
=======================================
CTSS clustering method: 
Number of TCs per sample:
=======================================
Consensus cluster information
=======================================
Number of consensus clusters:
Consensus cluster chromosome:
Consensus cluster start:
Consensus cluster end:
Consensus cluster strand:
Normalized tpm:
=======================================
Expression profiling
=======================================
Expression clustering method: 
Expression clusters for consensus clusters:
=======================================
Promoter shifting
=======================================
GroupX:
GroupY:
Shifting scores:
KS p-values (FDR adjusted):

Here, it is good to double check the reference genome and the sample labels. Also this gives us all the slots of the S4 object. So at this stage we have CTSS information as this is already the form of the data of the package. Only the tag count for each TSS for both samples, as the slot of Normalized tpm is empty. As well as the rest of the object.

Functions of CAGEr

With the following few standard CAGEr functions we will “fill in” these slots. And when full, show you how we can retrieve and export of these data.

Correlation plots

First, we will look at the data globally between the samples. CAGEr has a function do this for you (shown below) to look at the correlation between the samples and plots the png in the working directory. We fill first source an adjusted version of this function to control the output by using source (Creating a new generic function for ‘plotCorrelation’).

corr.m <- plotCorrelation(myCAGEset, samples = "all", method = "pearson") 
corr.m # correlation

            zf_512cells  zf_prim6
zf_512cells   1.0000000 0.4510721
zf_prim6      0.4510721 1.0000000

correlation

Normalization

Next, we need to normalize our data to adjust for example different library sizes to make them comparable. Here, we’ll use the power-law based normalization (see @balwierz).

librarySizes(myCAGEset)

[1] 4155747 7668701

# normalisation first plot to assess alpha
plotReverseCumulatives(myCAGEset, fitInRange = c(5, 1000), onePlot = TRUE) # plot
# insert alpha in normalize function:
normalizeTagCount(myCAGEset, method = "powerLaw",fitInRange = c(5, 1000), alpha = 1.20, T = 1*10^6)

Now the CTSS information the normalized tpm is added:

myCAGEset

Tag clusters of CTSSs

Let’s create tag clusters of CTSSs with the clusterCTSS function of CAGEr. CTSS in close proximity of each other give rise to functionally simialr set of transcripts within the same promoter elements. Tag clusters (TCs) are basically larger transcriptional units that correspond to individual promoters.
Here, we will use a simple distance measurement that individual CTSS can not be more than 20 bp apart. We also set the threshold at 1 which means that each TSS should have 2 or more tag counts in all the samples prior to clustering. Additionally, to remove low fidelity TSSs: no TC are called with a single TC if the normalized signal is below 5.

# Clustering of CTSS :)
clusterCTSS(object = myCAGEset, threshold = 1, thresholdIsTpm = TRUE, 
            nrPassThreshold = 1, method = "distclu", maxDist = 20, 
            removeSingletons = TRUE, keepSingletonsAbove = 5)

Keep in mind that these are determined per sample

Have a look again at the myCAGEset! The Tag cluster information slot is filled in.

Promoter width

Another feature we can check is the promoter width. That is, the width of each tag cluster per sample. CAGEr has a function that calculates the TC width by caculating cumulative distribution of tag signal per TC. Subsequently, promoter (TC) width is defined as the distance (in bp) between two quantiles generally set at 0.1 and 0.9 to capture 80% of tags.

# calculate cumlative CTSS distribution:
cumulativeCTSSdistribution(myCAGEset, clusters = "tagClusters")
# determine quantile positions:
quantilePositions(myCAGEset, clusters = "tagClusters", qLow = 0.1, qUp = 0.9)

We can then plot a histogram of the interquantile widths per sample with CAGEr too.

# we can plot the interquantile width with CAGEr too
plotInterquantileWidth(myCAGEset, clusters = "tagClusters",tpmThreshold = 3, qLow = 0.1, qUp = 0.9)

If you look at these two plots, you can see that the 512 cell stage more narrower TCs than prim6 which is what we would expect.

Consensus Clusters

CAGEr vignette:
Tag clusters are often sample-specific, thus can be present in one sample but absent in another. In addition, in many cases tag clusters do not coincide perfectly within the same promoter region, or there might be two clusters in one sample and only one larger in the other. To be able to compare genome-wide transcriptional activity across samples and to perform expression profiling, a single set of consensus clusters needs to be created. This is done using aggregate-TagClusters function, which aggregates tag clusters from all samples into a single set of non-overlapping consensus clusters:

aggregateTagClusters(myCAGEset, tpmThreshold = 5, qLow = 0.1, qUp = 0.9, maxDist = 100)

Expression - Self organising maps

As CAGE signals is essentially transcription, levels can be compared between groups. This can be done per CTSS, TC, or consensus clusters. CAGEr offers two methods to cluster gene expression, k-means and self-organising maps (SOM). Both require the number of expression cluster to be known in a priori.
Here we will perform the SOM algorithm at consensus cluster level. We set the threshold to at least 10 tpm (normalized) in at least one sample.

getExpressionProfiles(myCAGEset, what = "consensusClusters", tpmThreshold = 10,
        nrPassThreshold = 1, method = "som", xDim = 3, yDim = 2)
# in a nice plot:
plotExpressionProfiles(myCAGEset, what = "consensusClusters")

Promoter shifting

Promoter shifting is a method described by Haberle et al in 2014. They’ve shown that the same promoter can be used differently in different samples. This method has been implemented in CAGEr. Shifting can be detected between two individual samples (or between two groups of samples that are merged per group)

For all promoters a shifting score is calculated based on the difference in the cumulative distribution of CAGE signal along that promoter in the two samples. In addition, a more general assessment of differential TSS usage is obtained by performing Kolmogorov-Smirnov test on the cumulative distributions of CAGE signal, as described below. Thus, prior to shifting score calculation and statistical testing, we have to calculate cumulative distribution along all consensus clusters:

cumulativeCTSSdistribution(myCAGEset, clusters = "consensusClusters")

Calculating cumulative sum of CAGE signal along clusters...
    -> zf_512cells
    -> zf_prim6

Determine the shifting score:

Values of shifting score are in range between -Inf and 1. Positive values can be interpreted as the proportion of transcription initiation in the sample with lower expression that is happening “outside” (either upstream or downstream) of the region used for transcription initiation in the other sample. In contrast, negative values indicate no physical separation, i.e. the region used for transcription initiation in the sample with lower expression is completely contained within the region used for transcription initiation in the other sample.

To assess any general change in the TSS usage within the promoter region, a twosample Kolmogorov-Smirnov (K-S) test on cumulative sums of CAGE signal along the consensus cluster is performed. K-S test is performed to assess whether the two underlying probability distributions differ. To obtain a P-value, sample sizes that generated the ECDFs are required, in addition to actual K-S statistics calculated from ECDFs. These are derived either from raw tag counts (when useTpmKS = FALSE), or from normalized tpm values (when useTpmKS = TRUE). P-values obtained from K-S tests are further corrected for multiple testing using Benjamini and Hochenberg (BH) method and for each P-value a corresponding false-discovery rate (FDR) is also reported.

Here we will set the threshold of shifting score at 0.6 like the original paper as well as setting the FDR at 0.01:

This returns genomic coordinates (of conensusclusters), shifting score, and P-value (FDR) of the promoters, as well as the value of CAGE signal and position of the dominant TSS in the two compared (groups of) samples.

colnames(shifting.promoters)[7:8] <- paste(samples[c(4,11)], ".pos", sep = "") # to remember later which was group x and group y
write.table(shifting.promoters, "../Data/intermediate/ShiftedPromoters_sc06_fdr001_2samples.txt", col.names = TRUE,
            row.names = FALSE, quote = FALSE, sep = "\t")

Save myCAGEset as an intermediate file for if you want to revisit it later on

save(myCAGEset, file ="../Data/intermediate/CAGEobject_twoSamples_PowNom_allSlots.RData")

Export browser tracks from CAGEr

Export tables

Data handlin in R

Annotation of transcripts

Where do the cage signal

Dinucleotide plots

Differential Gene Expression

CAGE data can also be used to assess expression of the cTSSs. Here, we will use the R-package DESeq2 [@love_moderated_2014]. Originated for RNA-seq data but can also handle similar data from other assaya types (such as CAGE data). The vignette and reference manual can be displayed by running the code below or can be found here DESeq2.

browseVignettes("DESeq2")

Summary and goals of this practical

CAGEr

• Prepare the right data format from CAGE data for DESEQ2
• Export the data from CAGEr

DESeq2

• Normalise the data
• Differential expression

Follow up

• Gene annotation
• Gene ontology

1 Exporting data from a CAGEset object

Creating a count table input for DESeq2

We have worked until now with the two samples. However, for a differential expression analysis you will need more samples (replicates and/or more of the same condition). To this end, we’ll include two more samples to the mix to follow the more standard work-flow and generate p-values. The two additional samples are again from the same R package (ZebrafishDevelopmentalCAGE). The analysis will be on early stage expression vs later stage expression: zf_64cells and zf_512cells vs zf_prim6 and zf_prim20.

DESeq2 accepts matrices of read counts as input and this is exactly what we will export from CAGEr. It expects count data in the form of a matrix of integer values. The value in the i-th row and the j-th column of the matrix tells how many reads can be assigned to each cTSS site i in sample j. The values in the matrix should be un-normalized counts of sequencing reads.

Consensus site (row)	Sample 1 tpm	Sample 2 tpm	Sample 3 tpm	Sample 4 tpm
1	0	1	44	60
2	4	10	6	9

!! “The DESeq2 model internally corrects for library size, so transformed or normalized values such as counts scaled by library size should not be used as input”

The code below was run to produce the consensus clusters for the four samples as this would take to long to do on the day. The one important (different step) is shown below in the code. Don’t run the code (today).

# packages
require(ZebrafishDevelopmentalCAGE)
require(CAGEr)
# load data
data(ZebrafishSamples)
as.character(ZebrafishSamples$sample)
myCAGEset <- importPublicData(source = "ZebrafishDevelopment", dataset = "ZebrafishCAGE", 
group = "development", sample = as.character(ZebrafishSamples$sample[c(3:4,11:12)] ) )

# CTSS tag count
ctss <- CTSStagCount(myCAGEset)

#
# To keep using the raw counts in all downstream steps, the normalizeTagCount function of CAGEr should be used with the method set to "none". Note that normalizeTagCount function has to be applied to CAGEset object before moving to next steps.
normalizeTagCount(myCAGEset, method = "none")
#


# Clustering of CTSS: low fidelity cTSSs are removed (each cluster with only one cTSS signal < 5).
clusterCTSS(object = myCAGEset, threshold = 1, thresholdIsTpm = TRUE, 
            nrPassThreshold = 1, method = "distclu", maxDist = 20, 
            removeSingletons = TRUE, keepSingletonsAbove = 5)

# cumulative distribution and quantile positions
cumulativeCTSSdistribution(myCAGEset, clusters = "tagClusters")
quantilePositions(myCAGEset, clusters = "tagClusters", qLow = 0.1, qUp = 0.9)

# aggregate the clusters across the samples:
aggregateTagClusters(myCAGEset, tpmThreshold = 5, qLow = 0.1, qUp = 0.9, maxDist = 100)
save(myCAGEset, file = "../Data/provided/AggregatedTagClus_0109_4Samples.RData")

To be able to compare transcriptional activity across these samples, consensus clusters will be used for downstream analysis. This is provided to you in a .RData file. The following steps involve creating the count data table for each consensus cluster, merge it with the coordinates of the consensus, and to write a file with the output in the intermediate directory within the Data directory.

# package
require(CAGEr)

# load the data produced by the code above
load("../Data/provided/AggregatedTagClus_0109_4Samples.RData")

# create count tpm matrix per consensus cluster for each sample
count.df <- data.frame(consensusClustersTpm(myCAGEset))

# and the consensus coordinates (same order)
consensus.info <- consensusClusters(myCAGEset)

# create identifiers to link back
consensus.info$cons_clus_id <- paste("cid_",1:nrow(consensus.info), sep = "")
rownames(count.df) <- consensus.info$cons_clus_id

# save the combined info and count tpm as intermediate files 
# the order is the same so we can easily use cbind
count.consensus.info <- cbind(consensus.info[,-1], count.df)
write.table(count.consensus.info, "../Data/intermediate/CountTable_Consensus_4Samples.txt", col.names = TRUE, row.names = FALSE, sep = "\t", quote = FALSE)

# remove unnecessary files 
rm(myCAGEset, samples, consensus.info)

2 DESeq2 Data Analysis

First, we will make a DESeqDataSet object from the count table and add the “formula” which is the design of the analysis downstream (a linear model: ~ condition). In this example we are only using the variable (early vs late), however, if you want to add covariates in the model (e.g. batch) these would be extra columns in info.df (see below) and given in the model like : ~ column name covariate + variable. So the variable will be the last in the model.

# DESeq2 expects a matrix of count table:
count.matrix <- as.matrix(count.df) # if from previously saved file: as.matrix(count.consensus.info[,7:10])
head(count.matrix)

# the condition for the analysis in data.frame:
samples <- colnames(count.matrix)
info.df <- data.frame(condition = factor(x = c("early", "early", "late", "late"), levels = c("early","late")), row.names = samples)

The reason for specifying the levels of the factor is because by default, R will handle factor levels based on alphabetical order. It is good practice to get into, especially in analyses where it matters which level you want to compare against. Thus identifying correctly up- or down-regulated genes in this case.
The first thing to do now is to create a DESeqDataSet object that stores the count matrix and design:

require(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = count.matrix, colData = info.df, design = ~ condition)

Let’s first do a simple heatmap to check how similar (or dissimilar) the four samples are on a genome-wide transcription level. For this we will use the transformed data as performed by DESeq2. It can do three different types of transformation and here we will use the function rlog, which stands for regularized log. It transforms the original count data to the log2 scale by fitting a model with a term for each sample and a prior distribution on the coefficients which is estimated from the data (See DESseq vignette for more info).

# package
require(RColorBrewer)
require(pheatmap)
# Extracting transformed values:
rld <- rlog(dds, blind=FALSE)
# heatmap
sampleDists <- dist(t(assay(rld)))
sampleDistMatrix <- as.matrix(sampleDists)

rownames(sampleDistMatrix) <- rld$condition
colnames(sampleDistMatrix) <- NULL

colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)

pheatmap(sampleDistMatrix, clustering_distance_rows=sampleDists,
        clustering_distance_cols=sampleDists,col=colors)

Indeed, the samples are clustering according to the developmental time frame. There are more things you can do and check such as principal component analysis that are also described in the vignette.

The standard differential expression analysis steps are wrapped into a single function, ` DESeq. Results tables are generated using the function results, which extracts a results table with log2 fold changes, pvalues and adjusted pvalues. The text, condition treated vs untreated, tells you that the estimates are of the logarithmic fold change log2 (treated/untreated).

# the analysis
dds <- DESeq(dds)
# results
res <- results(dds)

Next, we want to know what are the lowest adjusted pvalues as these are of interest, which genes they represent, and gene ontology for patterns in the data.
First, we are adding the genomic coordinates of the cTSS still stored in:

count.consensus.info

and then reorder the results according to lowest p-value:

res.info <- cbind(count.consensus.info[,c(1:4,6)], data.frame(res@listData))
res.info <- res.info[order(res.info$padj),]
head(res.info)

An easy summary:

summary(res)  # orignal S4 output of DESeq

The amount cTSS differentially expressed (padj < 0.05):

sum(res$padj < 0.05, na.rm=TRUE)

Save the results in a table in the Intermediate directory:

# Dataframe for downstream
# save as intermediate file
write.table(res.info, "../Data/intermediate/DiffExpression_Consensus_4Samples.txt", col.names = TRUE, row.names = FALSE, quote = FALSE, sep = "\t")
# remove previous files
rm(res)

3 Gene annotation and Gene ontology

Let’s start by having a txdb agaub for our samples. The txdb was created with the code below and can be found in the provided data directory.

# dan rerio v7
require(GenomicFeatures)
require(AnnotationDbi)
txdb <- makeTxDbFromUCSC("danRer7", "ensGene")
saveDb(txdb, file = "../Data/provided/txdb_DanRer7.sqlite")

Load in the the txdb:

require(GenomicFeatures)
require(AnnotationDbi)
txdb <- loadDb("../Data/provided/txdb_DanRer7.sqlite")

Let’s run through a few things first and build a function for later purposes. First defining gene features to later assess in our differentially expressed cTSSs: promoters, upstream seq (5kb of promoter), exons, introns, and gene as annotated in ensemble. These will be in GRanges objects.

# Promoters (500bp window around refgene TSS)
promoters = trim(promoters(txdb, upstream=500, downstream=500))
# 5 kb upstream - 500bp 
upstream = trim(flank(promoters, 5000))
# exons
exons <- exons(txdb) 
# exons grouped by gene:
exons_gene <- reduce(exonsBy(txdb,"gene"))

# introns
introns = intronsByTranscript(txdb)
# genes
gene = genes(txdb)

To use the useful features of GenomicRanges we will have to convert our differential expression result to a GRanges object too:

# packages
require(BSgenome.Drerio.UCSC.danRer7)

total = res.info 
gr <- GRanges(seqnames = total$chr,
                ranges = IRanges(start = total$start,end = total$end),
                strand = total$strand,
                cons_clus_id = total$clus.gene, # the identifier
                seqlengths = seqlengths(Drerio))

References